function [model,momentStru,specStru,names]=...
    modelAnalysisSummaryAJ(modelHandle,param,model,sampleStru,priorStru,...
    reportStru,location,names,posStru)
% =========================================================================
% modelHandle:  handle to a function 
%
% param:        parameter vector 
%
% model:        structure 
%
%       model.solveopt 
%       model.addsol 
%       model.data 
%       model.trainvec 
%
% sampleStru    structure 
%
% *priorStru* structure with the prior 
%   .dist 
%   .mean 
%   .std 
%   .lbnd (lower bound) 
%   .ubnd (upper bound) 
%
% *reportStru* 
% 
% counterStru   innovations  
%               aZero 
% =========================================================================
cucd = pwd;
%% compareStru.pathName
% if ~exist('compareStru.pathName')
%     reportStru.pathName=['compare ',strrep(datestr(now),':','.')];
% else
%     reportStru.pathName=strcat(char(reportStru.pathName),' ',strdate);
% end
% %location.savePath=cr_dir(location.loadPath,reportStru.pathName); 
%% Load Priors
%loadPriorSummary
%% priorStru 
% priorStru.dist=prpar_dist; 
% priorStru.mean=prpar_lcone; 
% priorStru.std=prpar_lctwo; 
% priorStru.lbnd =prpar_lb; 
% priorStru.ubnd =prpar_ub; 


%% 1. reportStru: Assign default values for the moments 
[~,reportStru]=ch_field(reportStru,'NSimVec',[100 size(model.data,1)]);
[~,reportStru]=ch_field(reportStru,'NRep',100);
[~,reportStru]=ch_field(reportStru,'NCorrel',4);
[~,reportStru]=ch_field(reportStru,'NIRF',20);
[~,reportStru]=ch_field(reportStru,'NVarDec',8);
[~,reportStru]=ch_field(reportStru,'flagUnitIRF',0);

%% 2. Solve model 
disp(' '); 
disp('_______________________________________________'); 
disp('Solving Model'); 
[model.GG,model.RR,model.CONS,model.eu,model.cholV,model.ZZ,~,model.steady,~...
    ,names.steady,names.states,names.shocks]...
    =feval(modelHandle,param,model.solveOptions,model.addsol);
[names.observables,~]=extractObsNames(model.ZZ(:,:,1),names.states);

%% 3.  Create Report
names.stateNames = names.states;
names.shockNames = names.shocks;
names.obsNames = names.observables;

%% Top 
namesTab.outputFile = ['Report', strdate];
namesTab.model =location.root;
namesTab.specification=location.specification;
namesTab.subfolder=location.subfolder;
namesTab.param=names.param;

model.std =nan(length(model.param),1);
model.std(posStru.param.est) =reportStru.paramStd; 

if isfield(reportStru,'description')==false || isempty(reportStru.description)==true
    reportStru.description=[];
end
fclose('all');
fid = createReportSummary(location.savePath,namesTab,priorStru,model,reportStru.description);
%% reportPositions: obtain positions of all objects to report in reportStru
sampleStru.sampleVec = sampleStru.vector;
[reportStru,reportFlag]=reportPositions(reportStru,names,sampleStru);

%% statesSpectrum 
if isfield(reportStru,'stateSpectrum')&&~isempty(reportStru.stateSpectrum)
    fprintf(fid,'\\pagebreak \n');
    fprintf(fid,'\\section{Spectrum Decomposition} \n');
    disp('Begin Spectrum'); 
    [specStru.cell,specStru.table,specStru.xlsCell]=estimaSpectrum(modelHandle,param,model.solveOptions,model.addsol,...
        reportStru.stateSpectrum,reportStru,1000,location.savePath);
    clear formatStru;
    %     if length(specStru.table.colLabels) > 8
    %         formatStru.landscape = 1;
    %     else
    formatStru.landscape = 0;
    %     end
    fclose('all');
    cd(location.savePath);
    fid=fopen([namesTab.outputFile,'.tex'],'a+');
    cd(cucd); 
    latexTableCellSummary(fid,num2cprec(specStru.table.dataMat(:,:,end) ,3 ),...
        specStru.table.rowLabels,specStru.table.colLabels,formatStru,{'Spectral Decomposition'},...
        'Series','Spectrum');
    clear specStru.table;
    cd(location.savePath);
    save workspace;
    cd(cucd);
    disp('Spectrum Done');
    fprintf(fid,'\\pagebreak \n');
else 
    specStru=[]; 
end


%note: sampleSS == sampleStru
%flagStru = structureMerge(reportFlag,flagStru,1);
%clear tempflag;
%% momentBase: Compute all moments 
momentStru=momentsDrawsSummary(model.GG(:,:,end),...
    model.CONS(:,end),...
    model.RR(:,:,end),....
    model.cholV(:,:,end),...
    model.ZZ(:,:,end),reportStru);
%% 
cd(location.savePath);
% fid = fopen([namesTab.outputFile,'.tex'],'a+');
fprintf(fid,'\\pagebreak \n');

cd(cucd); 
%% Asymptotic Variance Decompositions: Observables
if reportStru.NVarDec > 0 
tableStru.dataMat=zeros([size(momentStru.varDecHorObs(:,:,1)) 4]); 
tableStru.dataMat(:,:,1)=momentStru.varDecHorObs(:,:,1); 
tableStru.dataMat(:,:,2)=momentStru.varDecHorObs(:,:,2); % should be 4
tableStru.dataMat(:,:,3)=momentStru.varDecHorObs(:,:,4); %should be 8
tableStru.dataMat(:,:,4)=momentStru.varDecObs; 

tableStru.rowLabels=names.obsNames; 
tableStru.colLabels=names.shocks.short; 
tableStru.corner={'Observables'};
tableStru.tableHeader={'Variance Decomposition 1 Period Horizon',...
    'Variance Decomposition 2 Period Horizon',...
    'Variance Decomposition 4 Period Horizon',...
    'Variance Decomposition Asymptotic'};
% latexStru.printPDF=false;
% % latexStru.fileName=['Variance Decomposition Observables',strdate];
% latexStru.fileName = namesTab.outputFile;
% latexStru.outPath=location.savePath;
% latexStru.fileExists=true;
% latexStru.maxCols=5;
% latexStru.maxRows=30;
% latexStru.orientation='portrait';
% latexStru.close = 0;
% latexStru.append = 1;
% latexStru.open = 0;

fprintf(fid,'\\section{Asymptotic Variance Decompositions} \n');
fprintf(fid,'\\subsection{Observables} \n');
% multiPageTable(tableStru,latexStru);

formatStru.landscape = 0;
% if length(tableStru.colLabels) > 8
%     formatStru.landscape=1;
% end

for counter=1:4
    latexTableCellSummary(fid,num2cprec( tableStru.dataMat(:,:,counter) ,3 ),...
        tableStru.rowLabels, tableStru.colLabels,formatStru,tableStru.tableHeader(counter),...
        'Series',['Decomposition F Horizon=',num2str(counter)]);
    %     if counter==2 || counter ==4
    %         fprintf(fid,'\\pagebreak \n');
    %     end
end
%% Asymptotic Variance Decomposition: States
ind = floor(log(reportStru.NVarDec)/log(2))+1;
tableStru.dataMat=zeros([size(momentStru.varDecHorStates(:,:,1)) ind]);
for i = 1:ind
    tableStru.dataMat(:,:,i) = momentStru.varDecHorStates(:,:,2^(i-1));
end
tableStru.dataMat(:,:,ind+1) = momentStru.varDecStates;
% tableStru.dataMat(:,:,1)=momentStru.varDecHorStates(:,:,1);
% tableStru.dataMat(:,:,2)=momentStru.varDecHorStates(:,:,2); % should be 4
% tableStru.dataMat(:,:,3)=momentStru.varDecHorStates(:,:,4); % should be 8
% tableStru.dataMat(:,:,4)=momentStru.varDecStates;

tableStru.rowLabels=reportStru.stateMom;
tableStru.colLabels=names.shocks.short;
tableStru.corner={'States'};
% latexStru.fileName=['Variance Decomposition States',strdate];
% latexStru.fileName = namesTab.outputFile;
% latexStru.maxCols =8;
% latexStru.maxRows =20;
% latexStru.orientation='portrait';
% latexStru.close = 0;
fprintf(fid,'\\pagebreak \n');
fprintf(fid,'\\subsection{States} \n');
% multiPageTable(tableStru,latexStru);
for counter=1:4
    latexTableCellSummary(fid,num2cprec( tableStru.dataMat(:,:,counter) ,3 ),...
        tableStru.rowLabels, tableStru.colLabels,formatStru,tableStru.tableHeader(counter),...
        'Series',['Decomposition F Horizon=',num2str(counter)]);
    %     if counter==2 || counter ==4
    %         fprintf(fid,'\\pagebreak \n');
    %     end
end

end
%% Impulse Responses
fprintf(fid,'\\pagebreak \n');
fprintf(fid,'\\section{Impulse Responses} \n');
names.rowNames = reportStru.stateIRF;
names.colNames = names.shocks.short;
handlesIRF = plotIRF(momentStru.IRFStates,[],names);
graphStru.handles = handlesIRF;
graphStru.names = names.colNames;
% graphStru.captions = names.colNames;
latexAddGraphsSummary(fid,graphStru,location.savePath); 
compileGraphsPDF(handlesIRF,'IRF',location.savePath); %removed one cd in code

%% Fit
fprintf(fid,'\\pagebreak \n');
fprintf(fid,'\\section{Fit} \n');
%% Volatilities
% Volatility simulated observables
in.modenames = {''};
in.varnames  = names.observables;
in.savename  = 'stdSimObs';
in.savepath  = location.savePath;
plot_stdhist(momentStru.STDSimObs, model.data,in);

% Volatility simulated states
if ~isempty(reportStru.stateMom)
    in.varnames = reportStru.stateMom;
    in.savename = 'stdSimStates';
    plot_stdhist(momentStru.STDSimStates, [],in);
end

% Asymptotic Volatility of Observations
in.modenames = {''};
in.varnames  = names.observables;
in.savename  = 'stdObs';
%in.savepath  = location.savePath;
plot_stdhist(momentStru.STDObs, model.data,in);

%% Asymptotic Correlations
% %Observables
[ ac_data, ~, ~] = accov(model.data,0:reportStru.NCorrel);
in.modenames = {'data','model'};
in.varnames = names.observables;
%in.savepath = outpath;
in.savename = 'crossCorrelationsObs';
gvec.handles=plot_all_autocorr_Summary(momentStru.corrSimObs,ac_data,in);
gvec.names=in.varnames;
%%
pdfToImport = {'stdSimObsDensity' 'stdSimStatesDensity' 'stdObsDensity' 'crossCorrelationsObs'};
for i = 1:length(pdfToImport)
    fprintf(fid,strcat('\\includepdf[page={-}]{','"',pdfToImport{i},'"','.pdf} \n'));
end
%% Close 
fprintf(fid,'\\end{document}');
fclose('all');



%% Counterfactual
% if nargin > 7 && isempty(counterStru)==false
%     disp('Performing Counterfactual');
%     if isfield(reportStru,'counter')==1 && ...
%             isempty(reportStru.counter.lessShocks)==false
%         disp('Deleting shocks');        
%         tempPos=cellposition( reportStru.counter.lessShocks, names.shocks.long );
%         counterStru.innovations(:,tempPos)=0;        
%     end
%     
%     %maybe this goes here
%     model.addsol.counter=str2func([func2str( model.addsol.funcKfilter ),'Counterfactual']); 
%     
%     [counterOut.yobs,statesCounter]=feval(model.addsol.counter,param,modelHandle,model.data,...
%         model.solveOptions,model.addsol,...
%         counterStru.innovations,counterStru.aZero);
%     counterOut.states=statesCounter(:,reportStru.stateDecompPos);
%     %counterOut.states=KFStru.smoothSt(:,reportStru.stateDecompPos);
% else
%     counterOut=[];
% end
% 
% 
% % Assuming we have priorStru AND model.stdLoopMat
% %% Smoother 
% if isempty(counterOut)==true && isempty(reportStru.stateDecomp)==false;
%     disp('Begin smoother');
%     [KFStru,logL]=feval(model.addsol.funcKfilter,param,modelHandle,model.data,model.trainVec,...
%         model.solveOptions,model.addsol);
%     KFStru.analysis.states=KFStru.smoothSt(:,reportStru.stateDecompPos);
%     KFStru.analysis.stateDecomp=KFStru.countSt(:,reportStru.stateDecompPos,:);
% else
%     KFStru=[];
% end
% 
% 
% 
